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Abstract 

We  have  developed  a  general  mathematical  model  to  represent  the  membrane-electrode  assembly  (MEA)  of  fuel-cell  systems.  The  model  is 
used  to  analyze  the  effects  of  various  polarization  resistances  on  cell  performance.  The  model  accommodates  arbitrary  gas  mixtures  on  the 
anode  and  cathode  sides  of  the  MEA.  Moreover,  it  accommodates  a  variety  of  porous  electrode  and  electrolyte  structures.  Concentration 
overpotentials  are  based  on  a  dusty-gas  representation  of  transport  through  porous  electrodes.  The  activation  overpotentials  are  represented 
using  the  Butler- Volmer  equation.  Although  the  model  is  general,  the  emphasis  in  this  paper  is  on  solid-oxide  fuel-cell  (SOFC)  systems  for  the 
direct  electrochemical  oxidation  (DECO)  of  hydrocarbons. 
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1.  Introduction 

Quantitative  models  can  be  valuable  in  the  interpretation 
of  experimental  observations  and  in  the  development  and 
optimization  of  systems.  The  objective  of  the  membrane- 
electrode  assembly  (MEA)  model  developed  herein  is  to 
provide  a  physically  based  model  that  can  be  used  to 
evaluate  the  effects  of  design  and  operating  alternatives. 
While  the  model  can  stand  alone,  the  software  is  written  in  a 
way  that  enables  incorporation  into  larger  system-level 
models  that  also  consider  fluid  flow  and  thermal  manage¬ 
ment  throughout  a  fuel-cell  system.  The  model  is  general  in 
the  sense  that  it  can  accommodate  a  variety  of  geometrical 
configurations.  It  can  also  represent  proton-conducting  or 
oxygen-ion-conducting  electrolytes.  The  discussion  and 
examples  in  this  paper,  however,  consider  solid-oxide 
fuel-cell  (SOFC)  system. 

The  great  advantage  of  SOFC  systems  for  highly  efficient 
electric  power  generation  lies  in  its  potential  for  direct  use  of 
hydrocarbon  fuels,  without  the  requirement  for  upstream 
fuel  preparation,  such  as  reforming  [1-4  .  A  typical  layout 
for  a  planar  design  of  the  SOFC  system  is  illustrated  in  Fig.  1 . 
The  membrane-electrode  assembly,  which  is  composed  of 
an  electrolyte  sandwiched  between  the  anode  and  the  cath¬ 
ode,  is  itself  sandwiched  between  metal  interconnect  struc- 
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tures.  The  fuel  and  air  channels  are  formed  in  the 
interconnect  structure.  The  electrolyte  is  a  dense  ceramic 
(e.g.  yttria- stabilized  zirconia  (YSZ)  or  gadolinia-doped 
ceria  (GDC)),  which  is  impermeable  to  gas  flow  but  is  an 

r\ 

oxygen-ion  (O  )  conductor.  The  composite  electrodes  are 
porous  metal-loaded  ceramics  (cermets)  (e.g.  Ni-YSZ  for 
the  anode  and  Sr-doped  FaMnCF  (FSM)  for  the  cathode), 
which  may  be  mixed  electronic  and  ionic  conductors  that 
can  extend  the  three-phase  boundaries  (TPB)  and  promote 
the  electrocatalytic  reactions. 

The  SOFC  system  involves  complex  transport,  chemical, 
and  electrochemical  processes,  with  its  operating  perfor¬ 
mance  strongly  affected  by  the  corresponding  transport 
resistances  and  the  activation  barriers.  Fig.  2  illustrates  some 
of  these  processes  at  the  MEA  level.  During  operation, 
oxygen  from  the  air  channel  transports  through  the  porous 
cathode,  whereupon  it  is  reduced  to  oxygen  ion  (O2-)  at 
gas-cathode-electrolyte  three-phase  boundaries.  The 
formed  oxygen  ion  at  the  cathode-electrolyte  interface  is 
then  transported  to  the  anode-electrolyte  interface  through 
the  ion-conducting  electrolyte.  At  the  same  time,  fuel  from 
the  fuel  channel  transports  through  the  porous  anode  to  the 
anode-electrolyte  interface,  where  it  is  oxidized  electroche- 
mically  at  the  gas-anode-electrolyte  three-phase  boundaries 
to  products.  The  products  (e.g.  H2O  and  CO2)  are  then 
transported  back  to  the  fuel  channel  through  the  porous 
anode  structure.  Transport  resistances  of  the  gas-phase 

r\ 

species  in  the  porous  electrodes  and  O  in  the  electrolyte, 
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Planar  fuel  cell  architecture  Porous  cermet  electrode  structure 


Fig.  1.  Typical  layout  of  a  planar  SOFC.  The  left-hand  panel  illustrates  the  system  at  the  scale  of  flow  channels  (in  mm),  while  the  right-hand  panel  illustrates 
the  small-scale  structure  (in  pm)  of  the  porous  cermet  electrode. 


as  well  as  the  activation  energy  barriers  for  electrochemical 
reactions  can  result  in  various  polarizations.  They  are  repre¬ 
sented  as  concentration  overpotentials,  activation  overpo¬ 
tentials,  and  ohmic  overpotentials.  The  main  contribution  to 
the  ohmic  polarization  is  from  the  transport  resistance  of 
02“  in  the  electrolyte.  The  concentration  polarization  is  due 
to  the  transport  resistance  of  the  gaseous  species  through 
porous  composite  electrodes,  and  the  activation  polarization 
is  related  to  the  charge-transfer  processes  at  the  electrode¬ 
electrolyte  interfaces.  These  various  polarizations  are  func¬ 
tions  of  both  operating  conditions  and  the  physical  proper¬ 
ties  of  the  cell.  The  operating  conditions  include 
temperature,  pressure,  and  fuel  and  oxidizer  concentrations. 
The  cell  properties  involve  materials,  macro-  and  micro¬ 
structures  of  the  electrolyte  and  the  composite  electrodes, 
which  include  porosity,  tortuosity,  permeability,  and  thick¬ 
ness  of  the  anode  and  the  cathode,  the  ionic  conductivity  and 
thickness  of  the  electrolyte,  the  active  area  and  activity  of  the 
electrode-electrolyte  interface. 

To  understand  these  various  resistances  and  enhance  the 
cell  performance,  parametric  studies  for  these  various  polar¬ 
izations  in  the  electrolyte  and  composite  electrodes  have 


Fig.  2.  Illustration  of  the  components  of  the  MEA  unit  cell. 


been  widely  made  through  experiments  [5,6],  theoretical 
analysis  [7],  and  numerical  models  [8-10].  The  theoretical 
analyses  and  numerical  models  for  the  activation  polariza¬ 
tions  require  an  understanding  the  elementary  thermal- 
electrochemical  reaction  mechanisms  and  the  microstruc¬ 
ture  of  the  electrode.  Due  to  the  complex  and  uncertain 
details  of  the  reaction  processes  (particularly  in  case  of 
hydrocarbons),  most  modeling  in  the  literature  has  focused 
on  using  H2  as  the  fuel.  Chan  et  al.  [8]  used  the  intrinsic 
charge-transfer  resistance  to  represent  the  activation  polar¬ 
ization.  Tanner  et  al.  [11]  have  derived  an  effective  charge- 
transfer  resistance  for  the  composite  electrode,  which  is  a 
function  of  the  microstructural  parameters  of  the  electrode, 
intrinsic  charge-transfer  resistance,  ionic  conductivity  of  the 
electrolyte,  and  the  electrode  thickness.  This  effective 
charge-transfer  resistance  model  was  later  used  to  analyze 
the  activation  overpotentials  [9,10].  Adler  et  al.  [12]  derived 
a  charge-transfer  resistance  model  for  the  composite  cath¬ 
ode,  which  is  a  function  of  the  tortuosity,  porosity,  interface 
area,  oxygen  self-diffusion  coefficient  and  oxygen  exchange 
coefficient.  The  activation  polarizations  for  a  particular 
electrode  may  also  be  measured  by  the  electrochemical 
impedance  spectroscopy  (EIS)  [5,6,13].  Previous  porous- 
media  gas-phase  transport  models  have  been  developed  to 
predict  concentration  overpotentials  [8-10].  These  include 
molecular  diffusion,  Knudsen  diffusion,  and  pressure-driven 
Darcy  flow.  However,  these  have  only  been  applied  in  fuel¬ 
cell  systems  with  binary  gas-phase  mixtures. 

The  model  developed  in  this  paper  represents  a  unit-cell 
structure  (Fig.  2),  considering  the  effects  of  various  over¬ 
potentials  on  the  cell  performance.  As  illustrated  in  Fig.  2 
the  cell  is  anode- supported,  with  a  very  thin  electrolyte  and 
cathode,  and  both  electrodes  are  porous  ceramic-metal 
composite  structures.  Channels  formed  in  an  interconnect 
structure  carry  the  flow  of  fuel  and  air.  Depending  on  the 
position  of  the  unit  cell  within  a  fuel-cell  system,  the 
composition  in  both  the  fuel  and  air  channels  varies  depend¬ 
ing  on  initial  composition  and  depletion. 
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2.  Model  formulation 


The  overall  performance  of  SOFC  systems  is  greatly 
affected  by  ohmic,  concentration,  and  activation  polariza¬ 
tions.  Thus,  the  operating  cell  potential  (£Ceii)  is  lower  than 
the  Nernst  potential  ( E ).  The  operating  cell  potential  (ECQ\\) 
can  be  formally  expressed  as: 


^cell  E  ^7conc,a  ^7act,a  h  ohm  ^conc,c  ^act,c 
^interface  h  leakage  ’ 


where  ^conc  a  and  f]conc  c  are  the  concentration  overpotentials 
at  the  anode  and  the  cathode  due  to  the  gas-phase  species 
diffusion  resistance,  ^act  a  and  f]act  c  the  corresponding  acti¬ 
vation  overpotentials,  f]ohm  the  ohmic  overpotential  in  the 
electrolyte,  77mterface  the  interface  overpotential  due  to  the 
contact  resistance  at  material  boundaries,  and  ^leakage  the  cell 
potential  loss  due  to  the  electric  current  through  the  electro¬ 
lyte.  Subsequent  sections  in  the  paper  develop  models  for 
each  of  these  overpotentials,  all  of  which  are  functions  of  the 
current  density  iQ. 


2.1.  Global  electrochemical  reaction 


A  global  electrochemical  reaction  can  be  written  as: 

v  ( Y  noXk\ + Vo'  ( Y  15  Y  vlak + Y  v'fXk’ 

\k=  1  J  \k=  1  /  k=  1  k=  1 

(2) 

where  Xk  is  the  chemical  symbol  for  the  kth  species  (which 
may  participate  as  a  fuel,  an  oxidizer,  or  both),  v /  and  vQr  the 
stoichiometric  coefficients  for  the  fuel  and  oxidizer  mix¬ 
tures,  and  v" k  and  v"  k  the  stoichiometric  coefficients  of 
the  kth  product  species  in  the  fuel  and  air  channels,  res¬ 
pectively.  Each  species  Xk  has  a  primary  identity  as  a  fuel, 
an  oxidizer,  a  product,  or  an  inert,  and  this  identity  plays 
an  essential  role  in  balancing  the  reaction  to  determine 
stoichiometric  coefficients  and  in  assigning  charge  trans¬ 
fer.  The  global  reaction  can  be  expressed  in  compact  form 
as: 

K  K 

YVk7k^YvkXk,  (3) 

k=  i  k=  i 

where  vk '  =  vf'nf,k  +  v0'n0,k  and  vk  =  v"k  +  v"k. 

Fuel  and  oxidizer  mixtures  are  separated  by  the  anode- 
electrolyte-cathode  ME  A  structure.  Thus,  in  writing  a  glo¬ 
bal  electrochemical  reaction  and  calculating  the  Nernst 
potential,  it  is  convenient  to  represent  the  fuel  mixture 
and  the  oxidizer  mixture  separately  as  reactants.  This 
separation  is  accomplished  by  writing  the  LHS  of  the 
reaction  as  two  summations.  The  fuel  channel  mixture 
composition  is  represented  in  terms  of  the  species  mole 
numbers  as  nf:k.  Similarly,  nQjk  is  used  to  represent  the 
mixture  of  the  composition  of  the  oxidizer  mixture  in  the 
air  channel  (e.g.  for  air,  nOQ0  =  0.21  and  nQ ^  =  0.79). 


The  RHS  of  the  global  reaction  is  also  separated  into  two 
terms.  Depending  on  whether  the  cell  is  an  oxygen-ion 
conductor  or  a  proton  conductor,  the  product  species  appear 
on  the  anode  or  cathode  side  of  the  ME  A.  Thus,  following 
reaction  balancing,  the  location  of  the  product  species 
affects  the  product  stoichiometric  coefficients  v" k  and 
v"*.  There  is  more  discussion  on  this  point  following  the 
discussion  on  reaction  balancing. 

Since  the  electrochemical  reactions  involve  charge  trans¬ 
fer  it  is  necessary  to  specify  the  charge  transfer  associated 
with  each  reactant  species.  For  this  purpose,  we  introduce 
two  sets  of  parameters  (zf,k  and  z0j),  which  specify  the 
charge  transferred  per  mole  of  the  kth  fuel  and  oxidizer 
species.  For  a  species  that  is  inert  relative  to  the  electro¬ 
chemical  charge-transfer  process,  z  =  0.  In  our  convention, 
the  product  species  (e.g.  H20  and  C02)  also  are  assigned 
z  =  0.  Consider  two  examples  for  assigning  z.  Since  02  in 

r\ 

the  air  channel  is  reduced  to  O  by  obtaining  four  electrons 

r\ 

(i.e.  02  +4e-  — >  20  ),  the  parameter  z0,o2  =  4.  Because 

r\ 

the  direct  electrochemical  reaction  of  CH4  with  O  at  the 
anode  to  produce  H20  and  C02  (i.e.  CH4  +  402~  — > 
C02  +  2H20  -j-  8e-)  releases  eight  electrons,  Zf,cH4  =  8. 
The  total  number  of  electrons  transferred  by  the  global 
electrochemical  reaction  is  represented  as: 

K  K 

He  =  Y  Vf 'nf,kZf,k  =  Y  V°'n°,kZo,k-  (4) 

k=  1  k=  1 

There  may  be  reasons  to  incorporate  oxidizer  species  into 
the  fuel  channel  or  fuel  species  into  the  oxidizer  channel.  For 
example,  oxygen  in  the  fuel  channel  may  promote  beneficial 
partial  oxidation  of  the  fuel  or  it  may  work  to  inhibit  deposit 
formation.  Such  species,  while  influencing  thermal  chem¬ 
istry  in  the  gas  or  on  catalytic  surfaces,  do  not  participate 
directly  in  the  electrochemical  reaction.  Thus,  they  are 
considered  inert  in  the  global  electrochemical  reaction. 
For  example,  if  02  (i.e.  an  oxidizer  species)  appears  as  part 
of  the  fuel  composition  (i.e.  in  the  first  term  of  the  global 
reaction  (2)),  then  Zf\o2  =  0.  The  added  oxygen  in  the  fuel 
channel  cannot  work  directly  to  produce  electric  energy — 
only  fuel  oxidization  with  oxidizer  from  the  cathode  channel 
can  produce  power.  Thus,  any  oxidizer  in  the  fuel  channel 
must  reduce  the  cell  potential  because  of  its  diluting  effect. 

Before  proceeding  to  evaluate  the  cell  potential,  the  global 
reaction  must  be  balanced  to  determine  the  stoichiometric 
coefficients.  Although  balancing  reactions  is  a  well  known 
process,  a  few  words  of  explanation  are  warranted.  A 
balance  equation  is  written  for  each  element  (e.g.  H,  C, 
O,  N),  requiring  that  the  reaction  conserves  the  element 
balance.  The  process  leads  to  a  system  of  linear  equations  for 
the  stoichiometric  coefficients  (vf',  v0r,  and  vk).  Because  of 
an  essential  linear  dependence  stemming  from  the  fact  that 
species  are  composed  of  elements  in  a  specified  way,  one  can 
arbitrarily  set  v/  =  1.  Then  depending  on  the  number  of 
product  species,  only  some  of  the  vk  are  non-trivial.  Spe¬ 
cifically,  the  number  of  product  species  (and  hence  the 
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number  vjj  coefficients)  is  one  fewer  than  the  number  of 
elements.  Once  the  overall  product  stoichiometric  coeffi¬ 
cients  v"  are  known,  then  either: 


or 


depending  on  the  type  of  cell  (i.e.  oxygen  ion  or  proton 
conducting)  and  hence  the  location  of  the  product  species. 

Inert  species  (i.e.  reactant  species  for  which  Zf,k  =  0  and 
Zoj  =  0)  are  initially  excluded  from  the  balancing  proce¬ 
dure.  Once  the  balance  is  completed,  the  inert  species  are 
reintroduced  with  the  same  stoichiometric  coefficients  on 
the  reactant  and  product  sides  of  the  global  reaction.  The 
entire  balancing  process  is  completely  automated  in  the 
software. 


2.2.  Nernst  potential  and  the  Nernst  equation 


Based  on  the  chemical  potential  balance  at  open  circuit, 
the  Nernst  equation  provides  that: 


E  =  E°  + 


RT 

n&F 


K 


k=  1 


Vo  V* 


where  E  is  the  Nernst  electric  potential,  E°  the  ideal  Nernst 
potential  at  standard  conditions  (po  =  1  atm),  pk  the  partial 
pressure  of  the  kth  species,  T  the  temperature,  R  =  8.314  J/ 
(mol  K)  is  the  universal  gas  constant,  and  F  =  96485.309  C/ 
mol  is  Faraday’s  constant.  The  ideal  Nernst  potential  at  the 
standard  conditions  is  given  as: 


A  G° 


where  A G°  is  the  change  in  standard-state  Gibbs  free  energy 
between  products  and  reactants  of  the  global  reaction  Eq.  (3). 
Specifically: 


AG° 


K 

5>f  ^44  "F  Vo  7lQ  k 

k=  1 


where  p°k  is  the  standard- state  chemical  potential  of  the  kth 
species.  The  standard-state  thermodynamic  properties  of 
ideal  gases  depend  on  temperature.  The  needed  thermody¬ 
namic  properties  are  readily  available  in  databases  such  as  in 
Chemkin  [14].  As  an  alternative  to  using  the  species  partial 
pressures,  the  Nernst  equation  can  also  be  represented  in 
terms  of  the  species  molar  concentrations  [Xk\  as: 


E  =  E°  +!E-Yfyi V  +  VoV*  -  V")  ln[Xt] 

Uet  k=l 


RT  ,  (RT 

^ - —  In  I  — 

neF  \po 


K 

y^pfikjc + voV* 

k=\ 


Fig.  3  illustrates  the  Nernst  potentials  for  three  fuels  as  a 
function  of  fuel  utilization  (at  T  =  750  °C  and  p  —  1  atm). 
Consider  the  situation  for  butane  where  the  global  reaction  is: 


C4H10  +  4  02  ->■  5H20  +  4C02. 


Fig.  3.  Nernst  potential  for  three  fuels-air  systems  as  a  function  of 
percentage  of  the  fuel  utilization.  As  the  fuel  is  “utilized”  it  is  converted  to 
stoichiometric  products  that  dilute  the  fuel  on  the  anode  side.  The  air  is  not 
depleted  in  this  system. 

If  the  fuel  is  50%  utilized  then  the  fuel  mixture  is  composed 
of  1/2  moles  of  C4H10,  5/2  moles  of  H2O,  and  4/2  moles  of 
CO2.  The  Nernst  potentials  are  highest  for  the  hydrocarbons, 
assuming  direct  electrochemical  oxidation.  The  figure  also 
shows  that  the  Nernst  potential  for  hydrogen  is  reduced  most 
as  a  function  of  fuel  depletion.  The  potentials  illustrated  in 
Fig.  3  are  calculated  presuming  that  the  oxidizer  is  pure  air 
(i.e.  no  depletion  of  the  oxygen  in  the  air). 

2.3.  Concentration  polarization 

For  open-circuit  conditions  (i.e.  zero  current  flow),  the 
species  concentrations  at  the  electrolyte  interface  (three- 
phase  boundary)  are  the  same  as  those  in  the  bulk  channel 
flow.  Denoting  the  molar  species  concentrations  in  the 
channels  as  [X^]*,  the  Nernst  potential  is  written  as: 

E*  =  E°  +^-J2(v(nu  +  v0'n0,k  -  v")  ln[Xt]* 

Het  k= 1 

H - t;  ln( — )  £(vfV  +  vo'n0,k  —  v" ) ■  (11) 

neP  \P0  J  k=\ 

However,  when  the  current  is  flowing  there  must  be  con¬ 
centration  gradients  across  the  electrode  structures.  This  is 
because  the  diffusion  processes  that  supply  the  species  flux 
to  the  electrochemical  reactions  are  driven  by  concentration 
gradients.  Consequently,  species  concentrations  at  the  three- 
phase  boundaries  [A7]s  are  different  from  the  bulk  concen¬ 
trations.  Evaluated  at  the  three-phase  boundaries,  the  corre¬ 
sponding  Nernst  equation  is  written  as: 

Es  =  E°  +ELj2{vf'n{tk  +  Vo  nok  -  v")  ln[X,]s 
Heti  k=  1 


(10) 
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The  potential  difference  associated  with  the  concentration 
variation  is  written  in  terms  of  a  concentration  overpotential 
as  follows: 


n  =  Ff  —  F: 

0  cone 


RT 

nPF 


K 


y>y+voV 


k= i 


The  equation  can  easily  be  rewritten  to  separate  the  anode 
and  cathode  contributions  as: 


K 


^1  cone 


RT^,  ,  ,M1 /W 

Z^(Vf  lni|xir 


nPF 


k=i 


+ 


}pF  *=i 


y^(vo'no,i  -  v"J  In 


M 


(14) 


Furthermore,  the  concentration  overpotentials  at  the  anode 
and  the  cathode  can  be  defined  separately  as: 


rj 


cone, a 


-  4,0  in( 


[Xk] 


s  /  5 


(15) 


i] 


conc,c 


k=  1  V 


(16) 


There  are  a  number  of  physical  processes  that  may  con¬ 
tribute  to  the  concentration  polarization.  These  include  gas 
species  molecular  transport  in  the  electrode  pores,  solution 
of  reactants  into  the  electrolyte  and  dissolution  of  products 
out  of  the  electrolyte,  and  the  diffusion  of  the  reactants/ 
products  through  the  electrolyte  to/from  the  electrochemical 
reaction  sites.  Typically,  however,  the  diffusive  transport  of 
the  reactants  and  products  between  the  flow  channels  and  the 
electrochemical  reaction  sites  across  the  electrode  is  the 
major  contributor  to  the  concentration  polarization.  For  an 
anode- supported  ME  A  structure,  the  concentration  polari¬ 
zation  at  the  cathode  is  usually  very  small  due  to  its  being 
very  thin  (e.g.  50  pm).  However,  the  concentration  polar¬ 
ization  at  the  anode  can  be  very  significant,  particularly  at 
the  high  current  density  and  high  fuel  utilization  rate. 


Kn  ^  1  bulk  diffusion  is  dominant  and  when  Knf^l 
Knudsen  diffusion  is  dominant.  When  Kn  ~  1,  which  is 
the  typical  situation  in  SOFC  electrodes,  both  bulk  diffusion 
and  Knudsen  diffusion  are  comparable  and  must  be  con¬ 
sidered  together.  The  dusty-gas  model  (DGM)  [15]  was 
developed  to  represent  multicomponent  transport  in  this 
situation. 

The  dusty-gas  model  can  be  used  to  describe  the  relation¬ 
ship  between  the  gas  species  concentrations  in  the  flow 
channel  and  concentrations  at  the  electrochemical  reaction 
sites  (at  the  electrode-electrolyte  interface).  The  DGM  was 
developed  initially  for  astrophysical  applications,  assuming 
that  the  pore  wall  consists  of  giant,  motionless,  and  uni¬ 
formly  distributed  particles  (i.e.  dust)  in  space.  Three  gas- 
phase  transport  mechanisms  are  considered:  molecular  dif¬ 
fusion,  Knudsen  diffusion,  and  viscous  flow.  Bulk  diffusion 
and  Knudsen  diffusion  are  combined  in  series  to  form  the 
total  diffusive  flux.  The  viscous  porous-media  flow  (Darcy 
flow)  acts  in  parallel  with  the  diffusive  flux.  There  are  three 
porous-media  parameters  that  appear  in  the  DGM.  They  are 
the  Knudsen  coefficient  (Kg),  the  permeability  constant  (Z?g), 
and  the  porosity-to-tortuosity  ratio  (</>g/ Tg),  all  which  can  be 
determined  experimentally.  The  DGM  can  be  written  as  an 
implicit  relationship  among  the  molar  concentrations,  molar 
fluxes,  concentrations  gradients,  and  the  pressure  gradient 
[16,17]: 


y-  [Xe]Nf  ~[Xk]Nf  | 
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In  this  relationship  Nf  is  the  molar  flux  of  species  k ,  [Xj  the 
molar  concentrations,  and  [XT]  =  p/RT  is  the  total  molar 
concentration.  The  mixture  viscosity  is  given  as  p  (which 
depends  on  the  mixture  composition  and  temperature)  and 
Dt,  and  D\  Kn  are  the  effective  molecular  binary  diffusion 
coefficients  and  Knudsen  diffusion  coefficients. 

The  effective  molecular  binary  diffusion  coefficients  in 
the  porous  media  D\(  are  related  to  the  ordinary  binary 
diffusion  coefficients  Dm  in  the  bulk  phase  [18]  as: 


2.3.1.  Dusty-gas  model 

The  concentration  polarization  is  caused  primarily  by  the 
transport  of  gaseous  species  through  porous  electrodes.  It  is 
affected  by  the  micro  structure  of  the  electrodes,  particularly, 
the  porosity,  permeability,  pore  size,  and  tortuosity  factor. 
There  are  four  major  physical  processes  that  affect  transport 
in  the  porous  media.  They  are  molecular  diffusion,  Knudsen 
diffusion,  surface  diffusion,  and  Darcy’s  viscous  flow.  The 
relative  importance  of  bulk  diffusion  (e.g.  Fickian)  and 
Knudsen  diffusion  depends  on  the  relative  frequency  of 
molecular-molecular  collisions  and  molecular-wall  colli¬ 
sions.  These  processes  can  be  characterized  by  the  Knudsen 
number  Kn  =  X / d,  where  X  is  the  mean  free  path  length  in 
the  gas  and  d  is  a  characteristic  pore  diameter.  When 


Dck(  =  ^DU.  (18) 

Tg 

Binary  diffusion  coefficients,  which  are  determined  from 
kinetic  theory,  may  be  evaluated  with  software  such  as 
Chemkin  [18].  Knudsen  diffusion,  which  occurs  due  to 
gas-wall  collisions,  becomes  dominant  when  the  mean  free 
path  of  the  molecular  species  is  much  larger  than  the  pore 
diameter.  The  effective  Knudsen  diffusion  coefficient  can  be 
expressed  as: 


(19) 


where  the  Knudsen  permeability  coefficient  Kg  =  rp</>g/ rg 
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(</>  is  the  porosity  and  Tg  the  tortuosity  factor),  Wk  are  the 
molecular  weights,  and  rp  the  average  pore  radius. 

Assuming  that  the  porous  electrode  is  formed  by  closely 
packed  spherical  particles  with  diameter  dp  (an  idealization), 
the  permeability  can  be  expressed  by  the  Kozeny-Carman 
relationship  as  [19]: 


3^2 


g  P 
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(20) 


Other  porous-media  situations  have  different  permeabilities. 

An  iterative  procedure  is  needed  to  determine  the  pressure 
and  concentration  gradients.  Assuming  that  the  species 
molar  fluxes  Nf  through  the  porous  electrode  are  known 
based  on  the  current  density  and  the  global  electrochemical 
reaction  at  the  electrode-electrolyte  interface,  the  pressure 
gradient  can  be  determined  as: 


— V/7  = 


EkNf  /Dyn 


(1  /RT)  +  {BJn)  E k[Xk]/D 
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This  expression  for  the  pressure  gradient  comes  from  a 
summation  of  Eq.  (17)  over  all  species  k,  recognizing  that 
summation  of  the  first  term  (representing  ordinary  diffu¬ 
sion)  vanishes  exactly.  The  concentrations  [Xk\  in  the 
denominator  are  evaluated  either  in  the  channel  or  at  the 
electrolyte  interface,  depending  on  the  direction  of  the  net 
molar  flux  Nk.  If  the  net  flux  is  toward  the  electrolyte  (as 
it  normally  is  in  the  cathode),  then  the  [Xk\  are  taken  as  the 
channel  concentrations.  However,  when  the  flux  is  toward 
the  channel  (as  it  usually  is  in  the  anode),  the  concentrations 
are  evaluated  at  the  electrolyte  interface.  The  iteration  is 
needed  because  the  electrolyte  interface  concentrations  are 
not  yet  known.  Assuming  an  initial  guess  for  the  concen¬ 
trations  [Xk\  (usually  the  channel  concentrations),  the  pres¬ 
sure  gradient  is  evaluated.  Then,  using  this  pressure 
gradient,  the  concentration  gradients  can  be  determined 
from  Eq.  (17)  as: 


-VM  =  £ 

£fk 


+ 


^k,Kn 


+ 


D t,Kn  h 


(22) 


The  concentrations  in  the  summation  on  the  RHS  are 
evaluated  as  discussed  above  depending  on  the  direction 
of  the  net  flux.  This  analysis  is  simplified  greatly  by  assum¬ 
ing  linear  profiles  of  pressure  and  concentrations  across  the 
electrode.  We  have  justified  this  assumption  by  solving  the 
two-dimensional  dusty-gas  model  over  a  large  range  of 
electrode  structures  and  fuel-cell  operating  conditions.  Once 
the  concentration  gradients  are  evaluated,  a  new  estimate  for 
the  concentrations  at  the  TPB  can  be  found  (assuming  linear 
concentration  profiles  and  a  known  electrode  thickness). 
With  the  new  estimated  concentrations  a  new  pressure 
gradient  can  be  found.  This  iteration  converges  very 
rapidly — typically  needing  only  two  or  three  iterations. 


Assuming  a  current  density  iQ  and  a  global  electroche¬ 
mical  reaction,  the  molar  flux  of  the  gas  species  through  the 
anode  and  cathode  can  be  evaluated  as: 


*ta  =  (Vf 
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(23) 

(24) 


Note  that  the  molar  flux  must  be  taken  as  zero  for  any  inert 
species  that  does  not  participate  in  the  electrochemical 
reactions  (e.g.  nitrogen  in  the  cathode).  Appendix  A  pro¬ 
vides  examples. 


2.4.  Activation  polarization  and  the  Butler-Volmer 
equation 


In  the  chemical  reaction  process,  the  reactants  must  usually 
overcome  an  energy  barrier,  that  is,  the  activation  energy.  In 
charge-transfer  electrochemical  reactions,  the  reactants  must 
overcome  not  only  a  thermal  energy  barrier  but  also  an  electric 
potential.  Thus,  a  portion  of  the  potential  cell  voltage  is  lost 
in  driving  the  electrochemical  reactions  that  transfer  the 
electrons  to  or  from  the  electrodes.  This  activation  polariza¬ 
tion  is  a  very  important  loss  when  the  electrochemical  reac¬ 
tions  are  controlled  by  the  slow  electrode  kinetics. 

Assuming  a  single  rate-controlling  reaction,  the  Butler- 
Volmer  equation  [20]  provides  a  relationship  between  the 
current  density  and  the  charge-transfer  overpotential  as: 


k) 


exp 


—  exp 


(25) 


This  equation  represents  the  net  anodic  and  cathodic  current 
due  to  an  electrochemical  reaction.  The  activation  over¬ 
potential  ^act  is  potential  difference  above  the  equilibrium 
electric  potential  between  the  electrode  and  the  electrolyte. 
The  factor  nf v  is  the  number  of  electrons  transferred  in  the 
single  elementary  rate-limiting  step  that  the  Butler-Volmer 
equation  represents.  It  is  quite  common  to  assume  n® v  =  1. 
The  anodic  and  cathodic  asymmetric  factors  (also  called 
charge-transfer  coefficients)  aa  and  ac  determine  the  relative 
variations  of  the  anodic  and  cathodic  branches  of  the  total 
current  with  an  applied  overpotential  ^act.  The  anodic  and 
cathodic  charge-transfer  coefficients  aa  and  ac  depend  on  the 
electrocatalytic  reaction  mechanism  and  typically  take 
values  between  zero  and  one.  Further,  they  are  constrained 
by  aa  +  ac  =  1.  The  exchange  current  density  io  is  the 
current  density  of  the  charge-transfer  reaction  at  the 
dynamic  equilibrium  electric  potential  difference  Eeq.  In 
other  words,  at  the  equilibrium  potential,  the  forward  and 
reverse  current  densities  are  equal  at  io.  A  high  exchange 
current  density  implies  the  reaction  proceeds  rapidly  upon 
varying  the  potential  from  its  equilibrium  value.  There  is  an 
activation  overpotential  for  the  anode  and  the  cathode. 
Consequently,  there  are  charge-transfer  coefficients  aa  and 
ac  for  both  the  anode  and  the  cathode. 
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There  are  two  potentially  interesting  limits  of  the  Butler- 
Volmer  equation.  At  very  high  surface  activation  overpo¬ 
tential  one  of  the  exponential  terms  in  Butler- Volmer  equa¬ 
tion  becomes  negligible.  The  resulting  Tafel  equation  is  a 
straight  line  on  a  semi-logarithmic  plot,  In (4)  ~  ^act.  Alter¬ 
natively  stated: 
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In  an  SOFC,  the  Tafel  equation  can  often  be  used  to  calculate 
the  activation  overpotential  at  the  cathode.  The  second 
limiting  case  is  realized  at  very  lower  activation  polariza¬ 
tion.  In  this  case,  the  following  linear  current  overpotential 
relationship  applies: 


RT  fiQ 

(aa  “1“  ac )neF  \Iq 

This  equation  can  often  be  applied  to  calculate  the  activation 
potential  at  the  anode  of  the  SOFC  system.  In  our  model 
neither  of  these  limiting  cases  is  actually  used  to  determine 
the  activation  overpotentials.  Rather,  the  Butler- Volmer 
equation  is  solved  iteratively  in  its  native  form. 

The  exchange  current  density  4  is  a  measure  of  the 
electrocatalytic  activity  of  the  electrode-electrolyte  inter¬ 
face  or  TPB  for  a  given  electrochemical  reaction.  In  fact, 
the  exchange  current  density  is  a  crucially  important  factor 
for  determining  the  activation  overpotential  f/act.  However, 
4  is  not  a  simple  constant  parameter,  but  its  value  may 
depend  on  the  operating  conditions  and  material  properties, 
including  concentrations  of  reactants  and  products  adjacent 
to  the  electrode,  temperature,  pressure,  the  microstructure 
and  electrocatalytic  activity  of  the  electrode,  and  even 
the  conductivity  of  the  electrolyte  [21].  Particularly  for 
composite  SOFC  electrodes  with  the  mixed  ionic  and 
electronic  conduction,  such  as  a  Ni-YSZ  for  the  anode 
and  a  LSM-YSZ  for  the  cathode,  the  electrochemical 
reaction  zone  can  be  extended  from  the  electrolyte-elec¬ 
trode  interface  into  the  electrode.  Thus,  the  charge-transfer 
resistance  can  be  reduced  and  the  exchange  current  density 
can  be  enhanced.  Furthermore,  the  exchange  current  den¬ 
sity  for  the  composite  electrode  becomes  a  function  of  the 
electrochemical  properties  of  the  electrolyte  and  electrode 
and  the  effective  TPB  length.  Thus,  micro  structural  proper¬ 
ties  of  the  composite  electrode,  such  as  the  grain  size  and 
volume  fraction  of  the  electrolyte,  will  affect  the  exchange 
current  density.  A  steady- state  current-voltage  experiment 
can  measure  4  at  a  particular  operating  condition.  How¬ 
ever,  it  is  very  difficult  to  obtain  a  general  analytical 


expression  for  4?  which  requires  detailed  understanding 
of  the  elementary  thermal  and  electrochemical  reaction 
mechanisms  occurring  at  the  electrode  and  the  microstruc¬ 
ture  of  the  electrode. 

In  our  model  we  provide  flexibility  for  the  exchange 
current  density  to  depend  on  the  reactant  concentrations 
as: 


K 

*o  =  *oIIV,  (29> 

k=  1 

where  is  a  constant  and  yk  and  Xk  the  reaction  order  and 
mole  fraction  (at  the  electrode-electrolyte  interface)  for  the 
kth  species,  respectively.  Although  this  formulation  provides 
flexibility  in  the  model,  it  is  still  far  from  a  general  descrip¬ 
tion  of  the  elementary  charge- transfer  reactions  steps. 
Furthermore,  it  is  difficult  to  determine  the  reaction  orders 
for  any  particular  system. 

At  the  cathode  of  an  SOFC  system,  the  gas-phase  oxygen 
reduction  and  ion  incorporation  into  the  electrolyte  is  a 
complex  process  that  involves  a  series  of  elementary  reac¬ 
tions  steps,  including  the  oxygen  adsorption  and  dissociation 
on  the  surface,  diffusion  of  the  intermediate  species  on 
the  surface,  charge  transfer,  and  oxygen  incorporation  into 
the  lattice  of  the  electrolyte.  Any  of  the  reaction  steps  may 
be  the  rate-limiting  reaction  for  the  overall  process,  depend¬ 
ing  on  the  operating  conditions  (temperature,  oxygen  partial 
pressure),  and  material  properties,  and  microstructure  of 
the  electrode  (catalytic  activity,  pore  size,  porosity,  tortuos¬ 
ity,  and  permeability).  For  the  LSM-YSZ  composite  cath¬ 
ode,  Kim  et  al.  [22]  have  indicated  that:  (1)  4  ~  Po?  ^ 
the  rate-determining  step  is  the  charging  step  of  adsorbed 
oxygen;  (2)  4  Pq{4  if  the  rate-determining  step  is  diffusion 
of  charged  surface  oxygen  to  the  TPB;  and  (3)  4  does  not 
depend  on  p0n  if  the  rate-determining  step  is  incorporation 
of  oxygen  ions  into  the  electrolyte.  For  the  Pt-YSZ  cathode, 
the  exchange  current  density  4  was  shown  to  be  proportional 
to  Pq{4  at  low  poo  and  p0]  4  at  high  p0l  [23].  Uchida  et  al. 
[21]  indicated  that  4  may  be  proportional  to  the  ionic 
conductivity  of  the  electrolyte  for  both  the  LSM-YSZ 
and  Pt-YSZ  cathodes  depending  on  the  operating  tempera¬ 
tures. 

The  situation  is  more  complex  at  the  anode  of  an  SOFC 
system.  Reaction  of  a  hydrocarbon  fuel  may  involve  direct 
electrochemical  oxidization  [24],  steam  and/or  dry  reform¬ 
ing,  pyrolysis,  and  deposit  formation.  The  general  thermal 
and  electrochemical  reaction  steps  may  include  dissocia¬ 
tive  adsorption  on  surfaces,  surface  and/or  charge-transfer 
reactions  among  the  adsorbed  intermediate  species,  surface 
and/or  bulk  diffusion  of  the  adsorbed  species,  and  deso¬ 
rption  of  products.  Even  for  hydrogen,  the  simplest  and 
most  widely  studied  fuel,  there  is  still  considerable  debate 
about  the  underlying  electrochemical  processes  [5,7,25- 
27].  The  rate-determining  reaction  step  of  the  hydrogen 
oxidation  at  the  anode  may  be  the  dissociative  adsorption 
of  hydrogen,  the  formation  of  hydroxyl,  a  charge-transfer 
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reaction,  or  the  desorption  of  water,  all  of  which  can 
depend  on  the  operating  conditions  (such  as  the  tempera¬ 
ture,  pressure,  and  species  concentrations),  macro-  and 
microstructures,  and  material  properties  of  the  electrode. 
Additionally  thermal  and  electrochemical  reactions  may 
occur  on  the  surface  and  in  the  bulk  of  the  metal  and/or 
ceramic  of  the  composite  cermet  electrode.  Thus,  there  is 
not  a  simple  means  to  represent  the  relationship  between 
the  exchange  current  density  and  the  species  concentra¬ 
tions  at  the  anode. 

There  are  ways  to  measure  the  exchange  current  density 
(or  infer  its  value  from  related  measurements).  Based  on  the 
Butler- Volmer  equation,  a  charge-transfer  resistance  Rct  can 
be  expressed  as: 


R~f  = 
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At  the  zero  activation  overpotential: 
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which  presents  a  relationship  between  the  exchange  current 
density  io  and  the  charge-transfer  resistance  Rct.  The  charge- 
transfer  resistance  can  be  determined  from  EIS  measure¬ 
ments  [6,13]. 


2.5.  Ohmic  polarization 

Sources  of  ohmic  losses  in  a  fuel  cell  are  the  resistance  to 
the  ion  flow  in  the  electrolyte  and  the  resistance  to  the 
electronic  flow  in  the  electrode.  Due  to  typically  high  metal 
loading  resistive  electric  conductivity  in  the  electrodes  is 
high,  which  limits  ohmic  losses  in  the  electrodes.  Thus, 
ohmic  polarization  in  an  SOFC  system  is  typically  domi¬ 
nated  by  ion  resistance  through  the  electrolyte.  Such  losses 
can  be  reduced  by  decreasing  the  electrolyte  thickness  and 
enhancing  its  ionic  conductivity. 

Since  electric  current  flow  in  the  electrolyte  and  electro¬ 
des  obey  Ohm’s  law,  the  ohmic  losses  can  be  expressed 
simply  as: 

^7 ohm  =  ^'e-^tot?  (32) 

where  Rtot  is  the  total  area-specific  cell  resistance,  including 
the  area- specific  resistances  in  the  electrode  Red  and  the  solid 
electrolyte  Re\.  The  ionic  conductivity  of  the  electrolyte  can 
usually  expressed  as: 

o-d  =  cro T~l  exp  (-  ,  (33) 

where  Ee\  is  the  activation  energy  for  the  ionic  transport. 
Thus,  the  specific  ohmic  resistance  of  the  electrolyte  with 


thickness  of  Le i  can  be  obtained  as: 

Ra  =  —  •  (34) 

O  e\ 

There  are  several  opportunities  for  interface  overpoten¬ 
tials  with  in  the  MEA  structure.  For  example,  anywhere 
two  materials  are  physically  bonded  together  there  can 
be  a  resistance  the  causes  an  interface  overpotential, 

^interface* 

If  the  electrolyte  has  electrical  conductivity,  there  is  a 
possibility  for  a  leakage  current  that  may  be  represented  as 
a  leakage  overpotential  ^ieaka2e.  In  our  model: 


^leakage 


(35) 


where  ^?eakage  is  the  leakage  overpotential  at  open  circuit  and 
/max  the  maximum  current  density  for  the  cell.  Once  the 
Nemst  potential  is  known,  ^{}eakage  can  be  observed  from 
measurement  of  the  open-circuit  potential.  The  value  of  /max 
can  be  determined  from  the  model  by  calculating  current 
density  at  which  the  cell  power  density  vanishes. 


3.  Solution  algorithm 

The  model  is  “solved”  by  determining  the  cell  voltage  for 
a  given  current  density,  or  vice  versa.  In  general,  the  cell 
potential  Eceii  is  expressed  as  the  difference  between  the 
Nemst  cell  potential  E  and  the  sum  of  all  the  relevant 
overpotentials,  which  depend  on  current  density  (Eq.  (1)). 
For  a  given  cell  configuration,  with  specified  geometry  and 
material  properties,  as  well  as  operating  conditions  (e.g. 
pressure,  temperature,  and  fuel  and  oxidizer  composition), 
each  of  the  overpotentials  must  be  evaluated  as  a  function  of 
current  density.  Since  the  overpotential  is  usually  an  implicit 
function  of  current  density  an  iteration  is  needed.  Most  of  the 
iterations  are  accomplished  using  a  Newton-Raphson 
method.  The  solution  procedure  follows  several  steps  as 
follows: 

•  Calculate  the  stoichiometric  coefficients  of  the  global 
electrochemical  reaction  by  solving  element  balance 
equations.  This  depends  on  the  specific  fuel  and  oxidizer 
composition. 

•  Based  on  the  stoichiometric  coefficients,  evaluate  the 
charge  transferred  nQ  and  the  Nernst  potential  E. 

•  Evaluate  the  species  molar  flux  through  the  anode  and  the 
cathode  (Nfik  and  N^k).  These  are  all  evaluated  explicitly 
from  the  current  density  and  the  stoichiometry  (Eqs.  (23) 
and  (24)). 

•  Calculate  the  species  compositions  at  the  electrode¬ 
electrolyte  interfaces  ([Xjs)  from  the  DGM  at  given 
species  molar  flux  and  species  compositions  in  the  chan¬ 
nels.  This  is  an  iterative  procedure. 

•  Evaluate  the  concentration  overpotentials  for  both  elec¬ 
trodes  (^conc  a  and  ^conc  c)  once  the  interface  concentra- 
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tions  are  known  (Eqs.  (15)  and  (16)).  This  is  an  explicit 
evaluation. 

•  With  a  specified  exchange  current  density  4  at  the  anode 
and  cathode,  invert  the  Butler- Volmer  equation  to  obtain 
the  activation  overpotentials  ^acta  and  rjactc.  This  is  an 
iterative  procedure. 

•  Calculate  the  ohmic  overpotential  rjohm  in  the  electrolyte 
(Eq.  (32)).  This  is  an  explicit  procedure. 

•  Determine  the  leakage  overpotential  ^ieakage  (Eq.  (35)). 
Assuming  that  ^i3eakage  is  specified,  the  value  of  zmax  must 
be  determined  iteratively.  The  maximum  current  density 
is  achieved  when  the  net  voltage  is  zero.  The  value  of  /max 
depends  on  the  cell  parameters  as  well  as  all  the  other 
overpotentials.  For  a  given  current  density  iQ,  evaluate  the 
leakage  overpotential. 

•  With  all  the  overpotentials  in  hand,  evaluate  the  cell 
voltage  Eceii . 

The  procedure  described  above  assumes  that  the  cell 
voltage  is  to  be  determined  once  a  current  density  is 
specified.  To  create  a  voltage-current  plot,  for  example, 
the  procedure  is  repeated  for  a  range  of  current  densities. 
When  the  MEA  model  is  incorporated  into  a  larger  system- 
level  fuel-cell  model,  it  is  usually  more  natural  to  specify  the 
cell  operating  voltage  and  determine  the  local  current  den¬ 
sity.  In  this  case  an  outer  iteration  is  implemented  around  the 
whole  procedure  outlined  above.  Despite  the  nested  itera¬ 
tions,  the  convergence  is  very  rapid. 

4.  Analysis  of  a  methane-fueled  SOFC 

The  objective  of  this  section  is  to  apply  the  model  to  a 
particular  anode- supported  SOFC  system  operating  with 
methane  as  the  fuel.  Table  1  lists  model  parameters  of  the 
unit  cell.  For  an  anode- supported  SOFC  system,  both  the 
electrolyte  and  the  cathode  are  much  thinner  than  the  anode. 
Thus,  it  may  be  anticipated  that  the  electrolyte,  ohmic  and 
cathode-concentration  overpotentials  will  be  relatively 
small.  The  transport  parameters  for  calculating  the  ohmic 
overpotentials  in  the  electrolyte  (Table  1)  are  based  on  the 
oxygen- vacancy  conductivity  [28]. 

Assuming  pure  CH4  as  a  fuel  and  air  as  the  oxidizer,  Fig.  4 
illustrates  the  predicted  cell  voltage  and  the  power  density  as 
a  function  of  the  current  density.  The  peak  power  is  seen 
to  be  about  0.4  W/cm2  at  about  0.4  V.  Fig.  4  also  shows 
the  contributions  of  all  the  overpotentials  as  shaded  regions. 
For  this  particular  system,  the  cathode  activation  represents 
the  largest  overpotential.  It  is  followed  by  the  anode-activa¬ 
tion,  ohmic,  and  anode-concentration  overpotentials.  The 
cathode-concentration  overpotential  is  so  small  that  it  does 
not  show  on  the  plot.  All  the  overpotentials  increase  with 
increasing  current  density.  The  sum  of  the  overpotentials 
reaches  the  Nemst  potential  (in  this  case,  about  1.15  V)  at 
the  highest  current  density  (4  ~  2.25  A/cm2)  where  the 
power  is  reduced  to  zero.  In  other  words,  the  cell  ceases 


Table  1 


Parameters  to  represent  a  CH4  SOFC  system 


Parameter 

Value 

Operating  temperature  ( T ,  °C) 

750 

Operating  pressure  ( p ,  atm) 

1 

Anode 

Thickness  (La,  pm) 

1000 

Porosity  (<fg) 

0.35 

Tortuosity  (ig) 

3.50 

Average  pore  radius  (rp,  pm) 

1 

Average  particle  diameter  (dp ,  pm) 

10 

Exchange  current  density  (i0,  A/cm2) 

0.40 

Charge-transfer  coefficient  (aa) 

0.50 

Number  of  electrons  (nf v) 

1 

Cathode 

Thickness  (Lc,  pm) 

50 

Porosity  (</>g) 

0.35 

Tortuosity  (ig) 

3.50 

Average  pore  radius  (rp,  pm) 

1 

Average  particle  diameter  (dp,  pm) 

10 

Exchange  current  density  (i0,  A/cm2) 

0.13 

Charge-transfer  coefficient  (aa) 

0.50 

Number  of  electrons  (nf v) 

1 

Electrolyte 

Thickness  (Le  1,  pm) 

20 

Activation  energy  of  O2-  (Ee\,  J/mol) 

8.0E+4 

Pre-factor  of  O2-  (cr0,  S/cm) 

3.6E+5 

Leakage  overpotential  (^max,  V) 

0.0 

to  produce  power  when  the  overpotentials  just  equal  the  cell 
potential. 

4.0.1.  Fuel  utilization 

Because  fuel  is  consumed  by  electrochemical  reaction  as 
it  flows  through  the  anode  channel,  the  fuel  stream  is  diluted 
by  product  species  (CO2  and  H20).  Furthermore,  as  fuel  is 


Fig.  4.  Cell  voltage,  overpotentials,  and  power  density  as  a  function  of  the 
current  density  for  CH4  as  a  fuel  in  an  SOFC  with  parameters  stated  in 
Table  1. 
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Fig.  5.  Cell  voltages  and  power  densities  as  a  function  of  the  current 
density  for  fuel  CH4  at  a  different  percentage  of  fuel  utilization.  The 
exchange  current  density  is  set  to  be  constant. 

used,  the  Nernst  potential  decreases  (Fig.  3).  The  various 
overpotentials  are  also  affected  by  the  fuel  depletion.  Fig.  5 
illustrates  the  cell  performance  for  situations  where  the  fuel 
is  50,  75,  and  85%  used.  In  this  example,  the  exchange 
current  densities  z'o  at  the  anode  and  cathode  are  set  to  be 
constants  (Table  1).  Furthermore,  the  cathode  channel  is 
presumed  to  be  pure  air  (i.e.  no  oxygen  depletion).  At  75  and 
85%  fuel  utilization  there  are  seen  to  be  critical  current 
densities  at  which  the  power  abruptly  terminates.  In  these 
cases  (approximately  1.57  A/cm2  for  the  75%  case  and 
approximately  0.88  A/cm2  for  the  85%  case),  the  anode- 
concentration  overpotential  becomes  so  large  that  no  fuel 
can  reach  the  anode  three-phase  boundary.  As  a  result  the 
cell  cannot  operate. 

The  cell  represented  by  Fig.  5  would  operate  well  at 
very  high  fuel  utilization.  For  example,  if  the  operating 
voltage  was  set  at  0.5  V  the  current  density  would  vary 
relatively  little  (approximately  between  0.6  and  0.8  A/ 
cm2)  with  corresponding  variations  in  power  density 
(approximately  between  0.42  and  0.32  W/cm2)  as  the 
composition  in  the  fuel  channel  varies  between  pure 
methane  and  a  highly  diluted  stream  where  the  fuel  is 
85%  depleted. 

Fig.  6  presents  the  result  of  another  simulation  where  all 
parameters  except  one  are  the  same  as  for  the  simulation 
shown  in  Fig.  5.  The  cell  represented  by  Fig.  6  presumes  that 
the  anode  exchange  current  density  is  first  order  in  the  fuel 
concentration.  That  is: 

k)  =  *o  Ach4  >  (36) 

where  XCh4  is  the  fuel  mole  fraction  and  /{]  =  0.65  A/cm2. 
As  the  fuel  is  depleted,  the  exchange  current  density  drops 
proportionally.  In  this  case,  the  results  are  strikingly  differ¬ 
ent.  As  the  fuel  is  depleted  the  current  density  and  the  power 
density  are  significantly  diminished.  The  anode- activation 
overpotential  dominates  the  losses.  It  would  likely  be 


Fig.  6.  Cell  voltages  and  power  densities  as  a  function  of  the  current 
density  for  CH4  at  a  different  percentage  of  fuel  utilization.  The  exchange 
current  density  is  assumed  to  be  the  first  order  in  the  fuel  concentration. 


impractical  to  operate  such  a  cell  at  a  fuel  utilization  greater 
than  about  50%. 

5.  Interpreting  experimental  observations 

Gorte  et  al.  [1,24]  have  measured  SOFC  performance  for  a 
variety  of  fuels  using  a  YSZ  electrolyte  and  a  Cu-ceria-YSZ 
cermet  anode.  Fig.  7  reproduces  some  of  the  reported 
measurements  together  with  results  from  our  model. 

Table  2  lists  parameters  used  in  our  MEA  model  to  fit 
the  experimental  data  for  the  case  of  H2  as  the  fuel.  Some  of 
the  physical  and  operating  parameters  are  reported  [1,24], 
but  we  have  estimated  others.  The  parameters  at  the  cathode 
are  chosen  such  that  the  concentration  overpotential  is 
very  small.  The  transport  parameters  for  calculating  the 


Fig.  7.  The  MEA  model  is  used  to  represent  a  range  of  experimental  data 
of  current-voltage  and  power  densities  for  an  SOFC  from  Gorte  et  al. 
[1,24]  by  varying  a  few  parameters. 
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Table  2 


Parameters  for  matching  the  H2  experimental  data  by  Gorte  et  al.  [1,24] 


Parameter 

Value 

Operating  temperature  ( T ,  °C) 

700 

Operating  pressure  ( p ,  atm) 

1 

Anode  (Cu-ceria-YSZ) 

Thickness  (La,  pm) 

400 

Porosity  (</>g) 

0.52 

Tortuosity  (ig) 

3.50 

Average  pore  radius  (rp,  pm) 

1.0 

Average  particle  diameter  (Jp,  pm) 

3.0 

Exchange  current  density  (ig,  A/cm2) 

0.0055 

Charge-transfer  coefficient  (aa) 

0.50 

Number  of  electrons  (ne) 

1 

Cathode  (LSM-YSZ) 

Thickness  (Lc,  pm) 

50 

Porosity  ( <f>g ) 

0.35 

Tortuosity  (ig) 

3.50 

Average  pore  radius  (rp,  pm) 

1.0 

Average  particle  diameter  (dp,  pm) 

3.0 

Exchange  current  density  (ig,  A/cm2) 

0.75 

Reaction  order  of  02 

0.5 

Charge-transfer  coefficient  (aa) 

0.5 

Number  of  electrons  (ne) 

1 

Electrolyte  (YSZ) 

Thickness  (Le  1,  pm) 

60 

Activation  energy  of  O2-  (Eeh  J/mol) 

8.0E+4 

Pre-factor  of  O2-  (cr0,  S/cm) 

3.6E+5 

Leakage  overpotential  (^max,  V) 

0.07 

ohmic  overpotentials  in  the  electrolyte  are  based  on  the 
oxygen- vacancy  conductivity  of  YSZ  [28].  To  match  the 
experimental  data  for  cases  of  CH4  and  C4H10,  the  para¬ 
meters:  porosity,  exchange  current  density,  reaction  orders 
and  the  leakage  overpotential  are  adjusted  as  indicated  in 
Table  3. 

As  illustrated  in  Fig.  7,  the  experimental  measurements 
can  be  represented  very  well  by  our  model  through  varying  a 
few  parameters.  We  must  hasten  to  caution,  however,  that  we 
claim  no  uniqueness  in  the  set  of  parameters  used.  Moreover, 
although  the  model  can  match  the  data,  it  is  difficult  to 
support  some  of  the  “best-fit”  parameters  on  physical 
grounds.  Nevertheless,  we  believe  that  one  can  gain  valuable 
insights  through  interpreting  experimental  observations  in 
the  context  of  a  quantitative,  physically  based  model. 

As  shown  in  Fig.  7,  the  measured  open-circuit  cell 
potential  for  the  H2  case  is  very  close  to  the  theoretical 
Nernst  potential.  However,  for  the  CH4  and  C4H10  cases  the 


measured  open-circuit  potentials  are  approximately  0.2  V 
lower  than  the  corresponding  theoretical  Nernst  potentials. 
Therefore,  we  have  assigned  a  leakage  overpotential  to 
capture  the  measured  open-circuit  potential  (Table  3). 

Close  inspection  reveals  that  the  current-voltage  curves 
for  these  three  fuels  have  substantially  different  functional 
forms.  The  varying  shapes  infer  that  the  overpotential  losses 
for  these  three  fuels  are  dominated  by  different  physico¬ 
chemical  processes.  For  the  H2  fuel,  the  nearly  straight-line 
shape  (for  the  entire  range  of  current  densities),  indicates 
that  the  major  polarization  may  be  due  to  the  ohmic  loss  at 
the  electrolyte  rather  than  the  activation  overpotential  and 
concentration  overpotential  at  the  anode.  For  CH4,  the  data 
shows  a  substantial  upward  curvature  at  low  current  density 
but  a  nearly  linear  function  at  high  current  density.  In  this 
case,  the  major  overpotential  loss  may  be  caused  by  the 
relatively  low  anode  activation.  In  the  case  of  C4H10,  the 
current-voltage  curve  is  again  functionally  different  from 
the  other  fuels,  particularly  at  the  high  current  density.  The 
relatively  lower  power  density  indicates  that  the  activation 
overpotential  is  very  important,  and  the  downward  curvature 
at  the  high  current  density  indicates  a  significant  anode- 
concentration  overpotential. 

Since  the  concentration  overpotentials  appear  to  be  rela¬ 
tively  small  for  cases  of  H2  and  CH4,  we  have  used  the 
measured  values  for  porosity  and  pore  size.  In  these  cases, 
we  have  varied  the  exchange  current  densities  and  reaction 
orders  to  fit  these  measurements.  However,  for  the  case  of 
C4H10  fuel,  the  anode-concentration  overpotential  appears 
to  dominate  the  potential  losses,  particularly  at  the  high 
current  density.  Consequently,  it  appears  that  the  gas-trans- 
port  resistance  through  the  anode  must  be  increased.  The 
most  direct  way  to  increase  transport  resistance  is  to  reduce 
the  anode  porosity.  A  very  low  porosity  of  <fi  ^  2%  is  needed 
to  reproduce  the  observed  strong  curvature  in  the  data.  Since 
the  measured  porosity  is  <fi  «  50%,  one  must  seriously 
question  the  physical  validity  of  the  porosity  needed  to  fit 
the  data.  It  is  more  likely  that  there  may  be  some  important 
physical  process  that  is  not  represented  in  the  model  as  it 
currently  stands.  For  example,  is  it  possible  that  there  is 
significant  catalytic  reaction  within  the  anode  structure?  If 
so,  perhaps  it  is  not  the  parent  butane  fuel  that  reacts 
electrochemically,  but  rather  some  reaction  products.  Cer¬ 
tainly  if  one  of  the  essential  underpinnings  of  the  model  (e.g. 
direct  electrochemical  oxidation  of  the  fuel)  is  incorrect, 
then  inferences  from  the  data-fitting  procedure  must  be 
questioned. 


Table  3 

Parameters  adjusted  to  match  the  experiment  data  by  Gorte  et  al.  [1,24] 

Fuel  type  Reaction  order  ( yk )  f  (A/cm2)  leakage  AO 

h2  ch4  c4h10  h2o  co2 


Hydrogen  (H2)  1.0  0.0  0.0  -1.0  0.0  0.520  0.0055  0.07 

Methane  (CH4)  0.0  1.0  0.0  -0.5  0.0  0.520  0.0045  0.20 

Butane  (C4H10)  0.0  0.0  1.0  -1.0  0.0  0.021  0.0090  0.23 
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6.  Summary  and  conclusions 

We  have  presented  a  physically  based  mathematical  model 
that  represents  a  unit  cell  of  a  fuel-cell  membrane-electrode 
assembly.  It  is  written  to  accommodate  general  descriptions 
of  fuel  and  oxidizer  mixtures.  It  considers  a  variety  of  anode 
and  cathode  overpotential  losses  that  include  activation, 
concentration,  ohmic,  and  leakage.  The  activation  overpo¬ 
tentials  are  described  in  terms  of  the  Butler- Volmer  equation. 
Concentration  overpotentials  in  the  porous  electrode  struc¬ 
tures  are  described  in  terms  of  a  dusty-gas  model.  Ohmic 
losses  may  stem  from  ion  resistance  in  the  electrolyte  or 
electrical  resistance  at  material  interfaces. 

The  model  can  be  exercised  to  produce  voltage-current 
relationships  for  a  certain  cell  configurations  and  specified 
operating  conditions.  These  simulations  are  valuable  in 
helping  to  understand  the  competing  physical  processes  that 
are  responsible  for  controlling  cell  performance.  Such 
understanding  can  assist  in  cell  design  and  optimization 
as  well  as  interpreting  experimental  observations. 

The  objective  of  this  paper  is  primarily  to  describe  the 
ME  A  model  and  the  underlying  assumptions.  The  model  can 
be  used  as  a  stand-alone  model.  However,  the  software  is 
written  to  be  incorporated  as  a  submodel  into  larger  system- 
level  models  that  also  consider  the  fluid  flow  and  thermal 
transport  for  a  full  fuel-cell  system. 

The  model  presented  here  assumes  direct  electrochemical 
oxidation  of  the  fuel  species.  Thus,  while  there  is  multi- 
component  transport  of  fuel,  oxidizer,  and  product  species 
within  the  electrode  structures,  there  is  no  homogeneous  or 
catalytic  reaction  among  these  species.  Although  there  is 
reported  experimental  evidence  of  DECO  [24]  for  hydro¬ 
carbon  fuels  in  certain  SOFC  systems,  there  is  also  reason  to 
anticipate  that  catalytic  reaction  could  occur  within  porous 
cermet  electrode  structures.  Thus,  in  future  work  we  intend 
to  extend  the  models  to  accommodate  this  possibility. 
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Appendix  A 


electrochemical  oxidation.  To  maintain  generality  there  is  a 
certain  level  of  abstract  nomenclature  that  must  be  used.  The 
objective  of  this  appendix  is  to  assist  understanding  and 
interpreting  the  nomenclature  by  working  through  specific 
examples. 

A.l.  C^HiQ-air  SOFC 


The  overall  electrochemical  reaction  is  stated  as: 


C4H10  + 


13 


2  x  0.21 


(0.21O2  +0.79N2) 


13  x  0.79 

4C02  +  5H20  +  - — — -N2, 


(A.l) 


2  x  0.21 

which  combines  the  electrochemical  reaction  in  the  anode: 


C4H10  +  130z 


4C02  +  5H20  +  26e_, 


(A.2) 


and  the  electrochemical  reaction  in  the  cathode: 

f  02  +  26e“  -►  1302".  (A.3) 

The  fuel  and  oxidizer  compositions  are  specified  by  setting 
the  non-zero  and  nQ^  as: 

ftf,c4Hi0  =  1,  (A.4) 

no,o2  =  0-21,  ft0,N2  =  0.79.  (A.5) 

Balancing  the  reaction  yields  the  non-zero  stoichiometric 
coefficients  as: 


Vf'  =  1, 

(A.6) 

13 

V°  2  x  0.21  ’ 

(A. 7) 

Vf,C02  =  4>  Vf,H20  =  5> 

(A. 8) 

„  13  x  0.79 

V°>N2  ““  2  x  0.21  ' 

(A. 9) 

The  non-zero  charge  transfer  per  mole  of  each  reactant 
species  are: 


£f,CH4  —  26,  z0,o2  —  4.  (A.  10) 

With  the  reaction  balanced  and  the  charge  transfers 
assigned,  the  number  of  electrons  transferred  by  the  global 
reaction  follows  from  Eq.  (4)  as: 

K 

ne  =  ff  vf'nitkZf,k  =  26.  (A.l  1) 

k=l 

The  concentration  overpotentials  are  expressed  as: 


This  paper  develops  and  discusses  a  general  formalism 
beginning  with  a  global  electrochemical  reaction  for  direct 


RT  1 

^conc,c  ^ 


\0 


\o 


(A.  13) 
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For  the  butane-air  case,  the  species  molar  fluxes  in  the  anode 
are: 


Vg  =  _ 

C4Hl0  26 F  ’ 


5 /e 

4F  ’ 


and  in  the  cathode  are: 


h_ 

4F  ’ 


4  4 

8F’ 

(A.  14) 


(A. 15) 


A.  2.  Fuel-mixture-air  SOFC 

This  example  considers  the  anode  channel  to  have  a 
mixture  of  fuel  (CH4),  products  (H20  and  C02),  and  oxygen 
(02).  At  some  particular  point  in  the  fuel-cell  system  the 
mole  ratios  in  the  anode  channel  are  presumed  to  be 
CH4:H20:C02:02=1. 0:0. 3:0. 15:0.05. 

This  oxygen  is  explicitly  added  to  the  anode  stream — it  is 
not  coming  through  the  MEA.  The  oxidizer  is  presumed  to 
be  air.  The  overall  electrochemical  reaction  can  be  repre¬ 
sented  as  follows: 


reaction  follows  from  Eq.  (4)  as: 

K 

IK-  =  vf'nf,kZf,k  =  8-  (A.24) 

k—\ 


The  concentration  polarizations  at  the  anode  and  cathode  are 
expressed  as: 


(A. 25) 


n 


conc,c 


(A. 26) 


The  species  molar  fluxes  in  the  anode  are: 


nch,  -  8F » 
nco2  =  -  gjF 


Ni  = 

U20 


4  F 

0, 


and  in  the  cathode  are: 


(A. 27) 


(CH4  +  0.3H2O  +  0.15CO2  +  0.0502) 

+  qjj  (O.2IO2  +  0.79N2) 

-►  (2.3H20  +  1.15C02  +  0.0502)  +  (5^T)0.79N2. 

(A.  16) 


Ne  =  — 

°2  4  f  ’ 


=  0. 


References 


(A. 28) 


The  half-cell  reactions  combine  the  electrochemical  reaction 
in  the  anode: 

CH4  +  402“  -►  C02  +  2H20  +  8e“ ,  (A.  1 7) 

and  the  electrochemical  reaction  in  the  cathode: 

202  +  8e“  — >  402-.  (A.  18) 

The  fuel  and  oxidizer  compositions  are  specified  by  setting 
the  non-zero  rif ^  and  n0^  as: 

^f,CH4  =  1,  ^f,H2G  =0.3,  ^f,co2  =0.15, 

«f,o2  =  0.05,  (A.  19) 

no,o2  =  0.21,  «0jn2  =  0.79.  (A. 20) 

Balancing  the  reaction  yields  the  non-zero  stoichiometric 
coefficients  as: 


V  =  1, 


v'  =  2 


0  0.21  ’ 


(A. 21) 


v 


// 

f,C02 


=  1.15,  vf,H20"  =  2.3,  vf,o,  =  0.05, 

2  x  0.79 


.// 


V 


.// 


o,N2 


0.21 


(A. 22) 


The  non-zero  charge  transfer  per  mole  of  each  reactant 
species  are: 


£f,CH4  —  8,  £0,o2  —  4.  (A. 23) 

With  the  reaction  balanced  and  the  charge  transfers 
assigned,  the  number  of  electrons  transferred  by  the  global 
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